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Введение. Отмечается, что использование в расчетах адаптивных сеток позволяет повысить точность и 
экономичность вычислительных алгоритмов, не увеличивая число узлов. Особенно эффективен такой подход 
при расчетах нестационарных задач. Цель данного исследования — разработка, построение и программная 
реализация методов построения расчетных двумерных оптимальных гранично-адаптивных сеток для областей 
сложной конфигурации с сохранением заданных особенностей формы и границы области. Применение таких 
методов способствует повышению точности, эффективности и экономичности вычислительных алгоритмов. 
Материалы и методы. Рассмотрена проблема автоматического построения оптимальной гранично-адаптивной 
сетки в односвязной области произвольной геометрии, топологически эквивалентной прямоугольнику. 
Получено решение для минимального набора входной информации: заданы граница области в физической 
плоскости и число точек на ней. Создание алгоритма и программы построения сетки базируется на модели 
динамики частиц. Это позволяет определять траектории движения отдельных частиц и исследовать динамику 
их парного взаимодействия в рассматриваемой системе. С помощью инструмента тазК отделяются внутренние 
и граничные узлы сетки, и это дает возможность определить скорости перемещения узлов с учетом специфики 
решаемой задачи. 

Результаты исследования. Разработанные методы построения оптимальной гранично-адаптивной сетки 
области сложной геометрии дают возможность решить проблему автоматического построения сетки в 
двумерных областях любой конфигурации. Для оценки результатов исследования алгоритма решена тестовая 
задача и визуализированы этапы решения. В виде рисунков показаны расчетная область тестовой задачи и 
работа функции расчета скорости перемещения внутренних узлов. Визуализация подтверждает преимущество 
такого метода построения сетки, при котором отделяются граничные и внутренние узлы. 

Обсуждение и заключения. Результаты теоретических и численных исследований важны как для изучения 
качественных свойств сеток, так и для развития методов построения расчетных сеток, позволяющих 


эффективно, с высокой точностью решать задачи численного моделирования. 
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Введение. Построение сетки изначально рассматривалось как необходимый вспомогательный шаг в 
решении других задач [1-3]. В частности, при решении задач вычислительной гидродинамики построение 
расчетной сетки — весьма трудоемкий и длительный процесс [4—6]. При исследовании акваторий реальных 
водоемов приходится иметь дело с областями с объективно предопределенной границей, проходящей вдоль 
береговой линии’. Граничные узлы прямоугольной сетки Декартовых координат, покрывающей водоем, могут 
точно не попадать на контур береговой линии. Следовательно, равномерная сетка должна быть очень мелкой, 
чтобы были приемлемыми вносимые ею погрешности задания береговой линии. Например, для Азовского моря 
двумерная сетка, как правило, содержит более полумиллиона узлов [7-11]. Отметим также, что при 
использовании прямоугольных сеток краевые условия ставятся в точках, смещенных от реальной границы, 


либо (при усечении граничных ячеек сетки) вблизи границы концентрируются неоднородности, связанные с 
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неравномерностью шагов по пространственным направлениям”. Учитывая сказанное, а также необходимость 
упрощения структур данных и алгоритмов их обработки, численное моделирование указанного типа задач 
целесообразно проводить на оптимальной гранично-адаптивной сетке [12—14]. 

В настоящей статье представлены алгоритм и программа, разработанные для построения 
четырехугольных оптимальных гранично-адаптивных расчетных сеток на основе метода динамики частиц в 
двумерной постановке. Этот мощный вычислительный метод позволяет представить узлы сетки как 
заряженные частицы и моделировать динамику системы, состоящей из огромного числа частиц (до миллиона). 
На примере тестовых задач демонстрируется эффективность работы алгоритма для областей со сложной 
границей. 

Материалы и методы. Описание метода построения 20 оптимальной гранично-адаптивной 
сетки. В Декартовой системе координат Оху введем область Ш). В области Ш) строится сетка 


«=|(х., у, = М, ] =Ьм, | по заданным координатам граничных узлов. Сеточная область разделяет 


область Р) на элементарные подобласти в виде четырехугольников. Узлы сетки перераспределяются вдоль 
координатных линий. Затем в нужных зонах, где требуется сгущение узлов, переменные х,у заменяются на 


сжимающие переменные ©, с помощью разделенных преобразований: 


(2,1) > (х,у): х=х(&,п), у=у(е,п), (1) 
(х,у) > (6,1): в=ё(%у), п=п(ху). (2) 
Таким образом, на плоскости (&,1) строится сетка @= {(& м, ), [= №, = ЬМ, | ‚ которая 


определяется функциями х = х(5,7), у = У(е,1). 


Построение расчетной сетки @ базируется на методе динамики частиц. Такая техника моделирования 
широко представлена в литературе, поэтому мы остановимся на ней лишь вкратце [15—17]. 

Узлы сетки @ представляем совокупностью частиц с зарядами 4, и массой т,‚, которые 
перемещаются в расчетной области О вдоль и вблизи ее границы. Частицы взаимодействуют друг с другом, 


причем силы взаимодействия имеют электрическую природу. 
На отдельно взятую 1-ю частицу со стороны /-й частицы по закону Кулона действует отталкивающая 


сила Р;. Абсолютная величина силы РЁ, определяется расстоянием между данными частицами, а ее вектор 


направлен противоположно радиус-вектору гу ‚ соединяющему 1-Й и /-Й заряды (рис. 1). 





Рис. 1. Схема взаимодействия частиц 
Траектории движения заряженных частиц определяют расположение узлов сетки. Обозначим 


координаты 1-й частицы т у, } ‚ а/-й — (2, |. Длина радиус-вектора гу ‚ определяющего перемещение узла 


(х, Е. ) к узлу Г 7, ) ‚ рассчитывается по формуле: 


в="(ж, ут, = (3) 
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Отметим особенность преобразования (1)—(2). Если расстояние г = 0, то узел Гы п, отталкивается от 


каждого из соседних. Они в свою очередь должны либо встать на предписанное расстояние, либо 
отодвинуться, притянувшись к соседнему узлу и освободив место для вновь вставляемого. 
Вид преобразования (1)—(2), сжимающего координаты х,у в зонах больших градиентов, определяется 


решением задачи. С этой целью в направлении осей Ох, Оу используется модельное уравнение, описывающее 


потенциал: 
| 
Е(х Ау; Уз й? ‚© я,)= ы (4) 
Вей 7’ ‚< т, 
где [ — коэффициент пропорциональности между потенциалом и расстоянием г в заданном узле, @ — 
некоторый параметр. 
Каждый узел стремится к снижению потенциальной энергии, а именно: 
УЕ (х, ухи» У } > т, =1, №, 7=1,М,, п=Ь М, т=Ь,М.. (5) 
Сила а У, ) связана с потенциалом Ё=А (х, С. „т следующим соотношением: 
^ о И 
Я (ху, = втаа (Е) = У (6) 


+2 
тп У 


При моделировании процесса взаимодействия подвижных частиц предполагаем, что узлы сетки, 
попавшие за пределы расчетной области, вынуждены переместиться в точку на границе области, расстояние до 
которой минимально. Схема взаимодействия между подвижной частицей и частицей на границе показана 


я 


на рис. 2. 






а) 6) 






Рис. 2. Схема взаимодействия частиц на границе расчетной области Ш): направление движения частиц внутри области 
О (а); направление движения частиц за пределами области /) (6) 


Технология численного интегрирования уравнений движения базируется на алгоритме, 
представленном в [17]. 

Описание программного средства для расчета сетки. Программа состоит из модулей, реализующих 
пять функций: управляющую, вызывающую расчетные функции заданное количество раз; три расчетные 
функции и функцию визуализации. Расчетные: функция, описывающая скорость перемещения узлов; функция 
перемещения узлов; функция проверки выхода узлов за границу области. 


Структура данных программы: А — входной массив размерами [№ ‚М | : 
С — массив визуализации; В,, В, — массивы размерами [и,т | ‚ описывающие расположение узлов; В, — 
массив размерами [и,т |, описывающий маски граничных условий; и, у — компоненты вектора скорости 


перемещения узлов; 1, / — счетчики; и, т — количество узлов по направлениям Ох, Оу соответственно; [ 


— коэффициент пропорциональности между потенциалом и расстоянием г в заданном узле; © — некоторый 
параметр (степень при г). Управляющая функция обнуляет массивы и задает начальное расположение узлов. 
Алгоритм, описывающий функцию скорости перемещения узлов. 
Входные массивы: Б., В,, ВБ, и параметр а (например, @ = 3). Выходные: и, у. 


1. Начало цикла по — переменным , ]. Значения счетчиков задаются равными 
1=0,...п—1, 7 =0,....т-1. 
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2. Отделение граничных и внутренних точек расчетной области. Если тазК = |1 — точка на границе 
области, переход на п. 3, если тазК = 2 — точка внутри области, переход на п. 4. 

3. Алгоритм, описывающий функцию скорости перемещения граничных узлов 

3.1. Обнуление массивов и, у. 


3.2. Начало цикла по переменным Й, Л. Значения счетчиков задаются равными 
Пе ныР Ва, Ее, В. 


3.3. Расчет расстояния от одного узла до другого: 





1 1 


о ) Е (Ву. — 


3.4. Проверка условия. Если и > 0, то переход к п. 3.5, иначе переход к п. 3.6. 


3.5. Расчет и, у для а=3: 


(Вх, —Вх, п) 
1,] и] 
ки Вт, 
„Л 1,] [03 11,] 
г 
(Ву. - Ву. я) 
О. 
1,] 1,] у° 1,1 


3.6. Наращивание счетчиков Й, /| и переход кп. 3.3. 
3.7. Наращивание счетчиков 1, / и переход к п. 3.1. 


4. Алгоритм, описывающий функцию скорости перемещения внутренних узлов 
4.1. Обнуление массивов и, у. 


4.2. Начало цикла по переменным Й, 1. Значения счетчиков задаются равными 
И =1-—1,...1+1 Л= 7-1... 7+1. 

4.3. Проверка условия. Если (Й-1)(Л-/)=0, то переход к п. 4.4, иначе переход к п. 4.8. 

4.4. Расчет расстояния от одного узла до другого: 
Ра ) + (Ву, —ВУ ). 
4.5. Проверка условия. Если выполняется и > 0, то переход к п. 4.6, иначе переход к п. 4.7. 
4.6. Расчет и, у для К = 0,005: 

и. <-и. —^ (Вх, — Вх, |} 






й, 7 
т зы м - (Ву, — Бул ). 
4.7. Наращивание счетчиков Й, /1 и переход к п. 4.4. 
4.8. Наращивание счетчиков 1, / и переход кп. 4.1. 


Алгоритм, описывающий функцию расчета перемещения узлов. 
Входные массивы: В,, В,, и, у и параметр / (например, / =30). Выходные массивы: В, В,. 


1. Начало цикла по переменным 17. Значения счетчиков задаются равными 
А У Вы И. 
2. Расчет массивов В,, В,: 
Вх. <Вх ТР 
Ву,, < Ву, +. 
3. Наращивание счетчиков 71, / и переход кп. 2. 


Алгоритм, описывающий функцию проверки выхода узлов за границу области. 
Входные массивы: В,, В,, А и параметр 4 (например, 4=3), описывающий размеры окна, в 


котором представлена расчетная область в случае выхода узла В,, В, за границу расчетной области. Выходные 


массивы: В, В,. 
1. Начало цикла по переменным 1,7. Значения счетчиков задаются равными 


о, БИО, 


2. Нахождение индексов в массиве А, соответствующих узлу (Вх.‚„Ву,, ) . 
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И <— |Вх, 





‚ Ло ВУ. 





3. Начальное значение расстояния до границы задается значением параметра Ч: 


г 24. 
4. Проверка условия выхода точки (п, Л) за границу расчетной области: если выполняется А, >0, 
иначе переход в п. 12. 
5. Начало цикла по переменным 12, 72. Значения счетчиков задаются равными 


12=-—а,...а, 12=-а,...а. 
6. Проверка точки (12+, УЕ Л) на принадлежность расчетной области: если выполняется 
нп.) л = О, Иначе переход в п. 10. 
7. Расстояние от узла (Вх, Ву.) до точки (12 +Й, 2+ Л) находится по формуле: 


2 


и] <— 





. . 2 . . 
Вх, , -2-й) + (Ву. -12-Л) 
8. Если 7[ > г ‚, выполняется п. 9, иначе переход к п. 10. 
9. Запоминаем точку расчетной области, ближайшую к узлу (Вх, Ву.) 
и<— 1, 13<«12+И, 1]3<72+Л. 
10. Наращивание счетчиков по переменным 12, /2 и возврат к п. 6. 
11. Смещение узла (Вх, Ву.) на границу расчетной области: 
Вх,, <—13, Ву, < ]3. 
12. Наращивание счетчиков 1, / и переход к п. 2. 


Алгоритм, описывающий функцию визуализации. 
Входные массивы: Б‚, В,, А. Выходной массив — С. 


|1. В массив визуализации помещаем массив АД. 


2. Начало цикла по переменным 1, /, К. Значения счетчиков задаются равными 


1, ] 


1=1,...п-3, ]=1,...т-2, К=|Вх |... Вх 





3. Наносим вертикальные линии: 


Ву. -ВУ, 
СЕ Ву + Зы. 


” Вх.,-Вх,, 
г 1,] 


(&- Вх, |< 255. 


4. Наращивание счетчиков Г, /, К и переход кп. 3. 


5. Начало цикла по переменным 1, 7, К. Значения счетчиков задаются равными 


=1,..,П-2, 1 =1,..,т-3, К= |Ву, 





..|Ву 


1, +1 





6. Наносим горизонтальные линии 


Ву. -Ву, 
вы Зы “.(к- Вх, | |< 255. 
Вх. ,-—Вх - 


И. 


СТА, Ву 


1] 


7. Наращивание счетчиков 1, /, К и переход кп. 6. 
Входной файл — ВМР. На нем черным цветом обозначена геометрия области, на которой строится 
сетка. Остальная область обозначена белым цветом. Исходный файл ВМР записывается в массив, при этом 


черный цвет соответствует значениям 0, а белый — 255. К выходной информации относятся массивы В,, В,, 


описывающие расположение узлов сетки, и массив С, в котором хранится геометрия исходной области с 
нанесенной сеткой. 
Результаты исследования. Результаты работы алгоритма демонстрирует решение тестовой задачи. 
Входные данные: исходная область вида, представленного на рис. З, а также расчетные данные 
п=12, т=14, 1=1,...,10, 7=1.....12, [=30, 4=3, а=3. 
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Рис. 3. Расчетная область тестовой задачи 


За начальное расположение узлов принимается расположение узлов сетки © без адаптации к границе 


исходной области. Визуализация работы функции перемещения внутренних узлов представлена на рис. 4. 


ОХУЕЕНЫЕНЕЕЯ.. 


1 
у 


ААА, 


х. { #+41 


2 


ЕЕ 





0 5 10 15 
Рис. 4. Работа функции расчета скорости перемещения внутренних узлов 


Результат построения четырехугольной гранично-адаптивной сетки @, покрывающей исходную 
область, получен на основе представленного алгоритма (рис. 5 а). Рис.56 показывает работу алгоритма 
программы для случая, когда не отделялись граничные и внутренние узлы. Явное преимущество сетки, 
представленной на рис. 5 а, состоит в том, что ее ячейки — выпуклые четырехугольники. Для сетки вида 5 6 
данное требование не выполняется. 





Рис. 5. Результаты работы алгоритма построения сетки (@) : изображение сетки @ , построенной при отделении граничных и 


внутренних узлов сетки (а); изображение сетки @ , построенной без отделения граничных и внутренних узлов сетки @ (6) 


Обсуждение и заключения. Предложена технология построения двумерных оптимальных гранично- 
адаптивных сеток, базирующаяся на методе динамики частиц. Разработан и опробован алгоритм численного 
расчета четырехугольных сеток для областей сложной конфигурации с сохранением заданных геометрических 
особенностей формы области и ее границы. На примере тестовой задачи исследовались возможности 
предлагаемого алгоритма. Одно из преимуществ данной программной реализации —Й возможность 
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автоматически расставлять узлы на границе расчетной области и получать выпуклые ячейки. Представленный 
подход обладает достаточной универсальностью и надежностью и может быть использован при триангуляции 
рассматриваемых областей. 
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